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ABSTRACT 

The relaxational behaviour of the bond-diluted two-dimensional Ising 
model below the percolation threshold is studied using Monte Carlo tech- 
niques. The non-equilibrium decay of the magnetization, M(t), and the 
relaxation of the equilibrium spin-spin autocorrelation function, C(t), are 
monitored. 

The behaviour of both C{t) and M{t) is found to satisfy the Kohlrausch 
law of a stretched exponential with the same temperature-dependent expo- 
nent. The Kohlrausch exponent does not appear to depend on the bond 
concentration. 

The results indicate that we are not yet in the asymptotic regime, even 
when C(t) and M(t) are less than 10 -4 . 
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1. INTRODUCTION 

The dynamics of disordered magnetic systems such as the dilute fer- 
romagnet have been studied recently both by analytic techniques [1-3] and 
Monte Carlo simulations [4-7]. 

In a ferromagnet [8], with a concentration p of sites or bonds, T c {p) 
denotes the critical temperature below which the system exhibits magnetic 
long-range order; this regime is indicated by F in figure 1 which depicts 
a schematic phase diagram for a diluted ferromagnet. Above T c (p) the 
system is a paramagnet (P). The boundary separating the two regimes 
intersects the zero-temperature axis at p c , the percolation threshold, i.e. 
T c (p = Pc) = 0. 

The highest possible critical temperature for ferromagnetic order in the 
disordered system permitted by a rare statistical fluctuation of the disorder 
over the entire system is called [9] the 'Griffiths temperature', Tq. For a 
diluted ferromagnet Tq is the transition temperature of the pure system, 
i.e. Tq — T c {p = 1). The regime T c {p) < T < Tq is usually referred to as 
the 'Griffiths phase' (G) [10]. 

The Hamiltonian of the dilute two-dimensional Ising ferromagnet is 
given by 

h = - JijS l s j (i) 

<i,j> 

where Si — ±1 are Ising spins, the summation runs over nearest neighbours 
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only and the exchange interactions are chosen from 

P(J ij ) = (l-p)6(J ij )+p6(J ij -l). (2) 

Setting both Boltzmann's constant and the exchange coupling to unity im- 
plies that T G = 2/ln(l + y/~{2)) for this system. 

A key quantity of interest in the study of the dynamics of the diluted 
ferromagnet is the equilibrium spin-spin autocorrelation function C(t), de- 
fined by 

C(t) = [< Si(t)Si(0) >], (3) 

where < . . . > and [. . .] denote thermal and disorder averages, respectively. 
The behaviour of C{t) depends on the phase under question. A novel form 
of dynamic scaling that involves t/^ 2+z results from an assumption [1] that 
the dynamics are dominated by compact clusters; here £ is the correlation 
length of order parameter fluctuations in the pure (p = 1) system and Z, 
the dynamical critical exponent, is also that of the non-dilute system. 

In this paper we present new data for the diluted two-dimensional Ising 
model below the percolation threshold. Our main objective is to try and es- 
tablish the validity of the predicted asymptotic behaviour of the equilibrium 
spin-spin autocorrelation function. 

In section 2 we give a brief review of the theory and the predicted 
asymptotic form of C(t). The simulation procedure is outlined in section 
3 and is followed by our results in section 4. The conclusion is given in 
section 5. 
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2. THEORY 

We review the theory which is primarily due to Bray [1,3]. In the 
following we shall concern ourselves with exponential factors only and ignore 
all power-law prefactors. Although the theory refers to site-dilution, the 
results are also applicable to bond-dilution with minor modifications. 

In Bray's [1] approach one writes 



where n a is the probability that a given site belongs to a cluster of type a 
(here a refers to size, shape and composition) and r a is the relaxation time 
for the cluster. Two distinct regimes are predicted. 

In the paramagnetic phase T >Tc, clusters of fully occupied sites (or 
bonds) dominate the long-time dynamics and a straight-forward variational 
argument leads to 



Note that for T = Tq it is expected that lnC(t) - -£ 2 /( 2 + z ) for all t since 
the correlation length (of the pure system) is then infinite. 

In the Griffiths phase (T c (p) < T < Tq) the behaviour depends on 
the form of the dynamics and it will be assumed that we are dealing with 
'model A' i.e. non-conserved dynamics [11]. C(t) is bounded below [2] by 
any subset of the terms forming the sum in equation (4). In particular, 




(4) 




for t <C i 2+z 
for t > i 2+z 



(5) 




(6) 
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where now / specifies just the cluster size. The probability that a given site 
belongs to a cluster of 'volume' I 2 is given to leading order by 

m ~ p l = exp[— I 2 ln(l/p)]. (7) 

The relaxation time, 77, is given by the Arrhenius result for the global 
reversal of the magnetization. Thus, 

n ~exp[al/T], (8) 

where a is the bulk surface tension at temperature T. It follows from 
equations (6) and (7) that 

C(t) > Y J ^v[-i 2 H^/v)-t/n\. (9) 
1 

The summation can be performed by substituting for 77 from equation (8) 
and then using the method of steepest descent for large t. As a result 

-lnC(t)/(ln£) 2 < A max for t » £ 2+z (10) 

with 

A max =(T/a) 2 \n{l/p). (11) 

Below the percolation threshold most clusters are ramified. The re- 
laxation time, r x , of a ramified cluster containing x sites is bounded above 
by that of a corresponding compact cluster with the same number of sites, 
that is 

t x < exp[ax 1/2 /T}. (12) 
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The number of x clusters per site below the percolation threshold and for 
large x is given by [12] 

n x ~ exp[— a(p)x], (13) 

where a(p) — > as p — > p c — . It follows from equations (4), (12) and (13) 
that an upper bound on C{t) is given by 

C{t) < exp[-a(p)x - te- axl/2/T }. (14) 

X 

The summation, which can once again be performed by the method of 
steepest descent for large t, leads to 

-lnC(£)/(ln£) 2 > A min for t > ^ 2+Z (15) 

with 

A min = {T/a) 2 a(p). (16) 
Thus, we can combine equations (10) and (15) and write 

lnC(t) - -A(lnt) 2 fort><e 2+z (17) 

where 

A-min — A < A maxi (18) 

and the amplitude A depends on the system parameters. 

A consequence of equation (17) is that the relaxation is non-exponential 
throughout the Griffiths phase. Hence, a dynamical phase transition occurs 
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at the Griffiths temperature from exponential to non-exponential relax- 
ation. Furthermore, the bounds on A, equation (18), confirm that this is 
the correct asymptotic form below the percolation threshold. As a(p), and 
hence A m i n , vanishes at p c , the upper bound is of no value for p > p c . It 
has, however, been argued that equation (17) holds throughout the entire 
Griffiths phase [1, 2, 13]. 

Computer simulations have been performed at the percolation thresh- 
old [4, 5, 7] and both slightly above [5, 6] and below [5] p c . The Monte Carlo 
data have been found to satisfy the Kohlrausch [14] law of a stretched ex- 
ponential which is given by 

lnC(t)~-Bt n 0<n<l (19) 

where B and n are temperature-dependent parameters. The numerical data 
can be fitted by the Kohlrausch form over a remarkable range of times, from 
very short times, up to the latest times achieved in the simulations. In 
particular, the stretched exponential form has been observed for the regime 

t^e +z 7]. 

Most of the simulation work so far has concentrated on the behaviour 
of the model at the percolation threshold. As argued above, equation (17) 
is believed to be the correct asymptotic form for the relaxation of C(t) for 
p < p c . This asymptotic form remains to be seen in the simulations for 
p > p c . Hence, it is important to establish that equation (17) does indeed 
describe the relaxation of C(t) for p < p c . There are very limited numerical 
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results available for systems below p c . The only other simulation that we 
are aware of is that of Colborne and Bray [5] who studied three-dimensional 
Heisenberg and Ising models slightly below p c . Good agreement with the 
theoretical predictions was found for the Heisenberg system. The data for 
the Ising system, on the other hand, were not at all conclusive. 

We, therefore, investigate the behaviour of a diluted Ising system below 
the percolation threshold. Our underlying aim is to look for the predicted 
asymptotic form of C(t). 
3. THE SIMULATION PROCEDURE 

The simulations were performed on a two-dimensional bond-diluted 
Ising model on a square 64 x 64 lattice. The Hamiltonian for our model is 
given by equation (1). For any given p and temperature, T, the nearest- 
neighbour interactions are chosen according to equation (2). We studied the 
system below the percolation threshold p c (= 1/2) for two bond concentra- 
tions: p = 0.1 and p = 0.25. We impose full periodic boundary conditions 
and perform simulations using conventional [15] Monte Carlo techniques. 

The correlation length, £(T,p), of the diluted system is given by [8] 

C\T,p) = ^ 1 (T, Pc ) + ^- 1 (0,p), (20) 

where £t and £ p are the thermal and percolation correlation lengths, re- 
spectively. Throughout our simulations we have [16] mm(£r , £p) <C the 
linear size of the lattice; the dominant clusters are not particularly large at 
the timescales of interest. Consequently, we do not expect our results to be 
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influenced by finite-size effects. 

At the beginning of the simulation all spins are aligned. The square 
lattice is regarded as two inter-penetrating sublattices and all the spins of a 
given sublattice are updated in parallel according to Metropolis dynamics 
with the transition probability given by 



where AE is the energy change involved. The thermal and disorder averages 
appearing in equation (3) are replaced by averages over time and samples, 
respectively. 

In the simulations we initially monitor the non-equilibrium decay of 
the magnetization, M(t), defined by 



and here N = 4096. The system is assumed to be in equilibrium once the 
magnetization has decayed to zero within the statistical error. We then 
measure the equilibrium spin-spin autocorrelation function from 



where t = refers to an equilibrium state. 

At asymptotic timescales the only mode which remains to be equili- 
brated is that associated with a global reversal of the magnetization. As a 



W(Si -> -Si) = min[l,exp(-AE/T)} 



(21) 




(22) 




(23) 
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consequence, Colborne and Bray [5] have argued that the non-equilibrium 
decay of M{t) should exhibit the same asymptotic behaviour as the relax- 
ation of the equilibrium C(t). The data for M(t) will be examined to see 
whether this is indeed the case. 

To obtain statistically reliable data we averaged over at least 10 5 sam- 
ples for each T and in some cases over 10 6 samples were used. 

4. RESULTS 

We now discuss our results. It should be noted that in our simulations 
we managed to obtain values of M{t) and C{t) an order of magnitude less 
than those achieved in previous work [4-7]. 

Our results are summarised in figures 2 to 6. Figure 2(a) gives a plot 
of our data for p = 0.1 for the spin-spin autocorrelation function in the 
form ln(-lnC(£)) against hit for both T > T G and T < T G . Clearly, for 
each temperature, there are two distinct linear regimes: an initial 'short- 
time' linear regime which becomes more prolonged as the temperature is 
lowered and a 'long-time' linear regime. In figure 2(a) we also display the 
lines of best fit which were obtained on assuming stretched exponential 
behaviour (equation (19)) for the 'long-time' regimes. The gradients of the 
plots decrease slowly with decreasing temperature. Hence, the Kohlrausch 
law is only satisfied if we assume a temperature-dependent exponent n c (T). 

Table 1 contains values of n c (T) extracted from figure 2(a) and these 
are plotted against T in figure 6 which also shows the values of other expo- 
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nents discussed below. 

Using the exact expression [17] for £ for T > Tq and Z = 2.076 ±0.005 
[18], it can readily be confirmed that t <C £ 2+z for T = 3.0 for the data 
displayed in figure 2(a). (For T = 4.0 we have t ~ £ 2+z for the displayed 
data.) This means that we satisfy the conditions under which the first 
part of equation (5) is expected to hold. It is interesting to note that 
whereas Bray [1] predicts n c = 2/(2 + Z) « 0.49, we find a temperature- 
dependence that ranges from n c (T = 3.0) = 0.723 ±0.020 to n c (T = 1.0) = 
0.513 ± 0.004. In our earlier work [4] on the diluted Ising model exactly 
at the percolation threshold we found the Kohlrausch exponent to range 
from n c (T = 3.0) = 0.704 ± 0.005 to n c (T = 1.5) = 0.487 ± 0.001, with 
n c (T = Tq) = 0.641 ±0.006. These values are comparable with those found 
in the present study below p c . 

In figure 2(b) we re-plot the data for T < Tq as ln(— lnC(t)) against 
ln(lnt). This time the exact expression [17] for £ for T < Tq implies that 
the requirement of equation (17) is met. The straight lines, which are guides 
to the eye, have gradient 2. It is clear from figure 2(b) that we have yet to 
attain the asymptotic regime. 

We now turn to our data for the decay of C(t) for p = 0.25. The 
data, which is displayed in figures 3(a) and 3(b), was analysed in a similar 
way to that for p = 0.1. Once again, the data in the 'long-time' regimes 
in figure 3(a) may be fitted to a stretched-exponential only if we assume a 
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temperature-dependent exponent n c (T); values of the exponents extracted 
from the fits in figure 3(a) are given in table 2 and displayed in figure 6. The 
Kohlrausch exponents would not appear to depend on the concentration, p. 

Figure 3(b) gives a plot of ln(— lnC(t)) against ln(lnt) of the data for 
T = 2.0. It is rather puzzling to note that the straight line of gradient 
2 would appear to indicate that we have already attained the asymptotic 
limit given by equation (17) even though our data do not yet meet the 
requirement of t ^> £ 2+z for this particular temperature! 

We now discuss our results for the non-equilibrium decay of the mag- 
netization. We plot our data for M(t) for p = 0.1 in figures 4(a) and 4(b) 
which should be compared with figures 2(a) and 2(b), respectively. We 
see that the plots are very similar. The stretched-exponential exponents, 
n m (T), extracted from the 'long-time' regimes of figure 4(a) are given in ta- 
ble 1 and displayed in figure 6. We see that the exponents n c {T) and n m (T) 
cannot be distinguished within the uncertainty of the data, especially for 
the higher temperatures. This is to be contrasted with the earlier results 
[5,7] which would appear to imply that n m (T) > n c (T) for the diluted Ising 
model exactly at the percolation threshold. 

The data for M{t) for p = 0.25 are presented in figures 5(a) and 5(b), 
respectively. The values of n m (T) extracted from figure 5(a) are given in 
table 2 and displayed in figure 6. Figure 5(b) is very similar to figure 3(b) 
in that the asymptotic behaviour predicted by equation (17) is seen even 
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though t <C £ for the data exhibited. 

We see from figure 6 that the value of the stretched-exponential expo- 
nent appears to depend only on the temperature. 

5. CONCLUSION 

We have presented new Monte Carlo data on the relaxation of the 
diluted two-dimensional Ising model below the percolation threshold. 

We find that the non-equilibrium decay of the magnetization and the 
relaxation of the equilibrium spin-spin autocorrelation function can both 
be described by the Kohlrausch law of a stretched exponential with a 
temperature-dependent exponent n(T). This result is very similar to previ- 
ous Monte Carlo work exactly at the percolation threshold [4-7]. However, 
unlike the earlier work [5-7] for the diluted Ising model exactly at p c , we 
also find evidence that the decay of M{t) and the relaxation of C(t) can be 
described by the same temperature-dependent Kohlrausch exponent. 

The exact asymptotic behaviour predicted by a recent theory [3] is not 
seen even though we managed to obtain values of M(t) and C(t) more than 
an order of magnitude less than in previous simulations. 

It may be that the predicted asymptotic decay occurs only when M(t) 
and C{t) are outside any observable regime. 
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Table 1 



T 


n c 


n m 


1.0 


0.513 ±0.004 


0.538 ±0.004 


1.5 


0.597 ±0.012 


0.649 ±0.018 


1.75 


0.641 ±0.021 


0.661 ±0.011 


1.90 


0.686 ±0.027 




3.0 


0.723 ± 0.020 




4.0 


0.779 ±0.017 


0.831 ±0.017 



Table 2 



T 


n c 




2.0 


0.670 ±0.013 


0.703 ±0.010 


2.27 


0.679 ±0.011 




3.0 


0.656 ±0.029 


0.785 ± 0.027 


4.0 


0.796 ±0.021 


0.827 ±0.029 
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TABLE CAPTIONS 



Table 1 

A table of the Kohlrausch exponents n c and n m for the decay of the equilib- 
rium spin-spin auto-correlation function and the non- equilibrium relaxation 
of the magnetization, respectively, for p = 0.1. 



Table 2 

A table of the Kohlrausch exponents n c and n m for the decay of the equilib- 
rium spin-spin auto-correlation function and the non-equilibrium relaxation 
of the magnetization, respectively, for p = 0.25. Note that Tq ~ 2.27. 
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FIGURE CAPTIONS 



Figure 1 

A schematic phase diagram for a diluted ferromagnet: paramagnetic (P), 
ferromagnetic (F) and Griffiths (G) phases are indicated. Tq is the Griffiths 
temperature. 

Figure 2(a) 

A plot of ln(— lnC(t)) against Int for p = 0.1. The straight lines are best 
fits assuming a stretched exponential form for the 'long-time' linear regimes 
(see text). Figure 6 displays the gradients n c (T) obtained from this plot. 
The temperatures shown are:o, T = 4.0; □, T = 3.0; o, T = 1.9; A, T = 1.75; 
m,T = 1.5; and n,T = 1.0. 

Figure 2(b) 

A re-plot of the data displayed in figure 2(a) as ln(— In C{t)) against ln(lnt). 
The straight lines, which are guides to the eye, have gradient 2. The tem- 
peratures shown are:o, T = 1.9; A, T = 1.75; ;T = 1.5; and □, T = 1.0. 

Figure 3(a) 

A plot of ln(— lnC(t)) against lnt for p = 0.25. The gradients n c (T) ex- 
tracted from the straight line fits are displayed in figure 6. The temperatures 
shown are: A, T = 4.0; o, T = 3.0; o, T = T G ; and □, T = 2.0. 
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Figure 3(b) 

A re-plot of the data for T = 2.0 shown in figure 3(a) as ln(— In C{t)) against 
ln(lnt). The straight line has gradient 2. Note that, despite the reasonable 
fit, we have yet to attain the limit t ^> £ 2+z at this temperature. 

Figure 4(a) 

A plot of ln(— lnM(t)) against Int for p = 0.1. The gradients n m (T) 
extracted from the straight line fits over the 'long-time' regimes are shown 
in figure 6. The temperatures displayed are: A,T = 4.0; o, T = 1.75; 
□,T = 1.5; o,T = 1.0. 

Figure 4(b) 

A re-plot of some of the data displayed in figure 4(a) as ln(— In C{t)) against 
ln(lnt). The straight lines are guides to the eye and have gradient 2. The 
temperatures shown are: o, T = 1.75; □, T = 1.5; o, T = 1.0. 

Figure 5(a) 

A plot of ln(— lnM(t)) against lnt for p = 0.25. The slopes of the straight 
lines, n m (T), are displayed in figure 6. The temperatures shown are: o, T = 
4.0; □,T = 3.0; o,T = 2.0. 

Figure 5(b) 

A re-plot of the data for T = 2.0 displayed in figure 5(a). The straight line 
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has gradient 2. Note that we are not yet in the limit t ^> £ + ; compare 
with figure 3(b). 

Figure 6 

A plot of the Kohlrausch exponents n CjTn (T,p) as extracted from the ear- 
lier figures against the temperature. Symbols: o, n c (T, 0.1); o, n m (T, 0.1); 
□, n c (T, 0.25); □, n m (T, 0.25). 
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